S7.3) Big data approximations with discrepancy

This notebook contains the code of the paper "Bayesian Calibration of Imperfect Computer Models using Physics-Informed Priors". The models are fitted in rstan and the code is available in the folder "STAN/Approximations".

Load packages

knitr::opts_chunk$set(
  message=FALSE,
  warning = FALSE,
  comment = '', 
  fig.width = 6, 
  fig.height = 4,
  fig.align = 'center'
)
getwd()

[1] "/Users/michais/Desktop/BCPI_codes_new"

Reality and modelling choice

\[\begin{align} \mathcal{R}: \quad & \frac{d P(t)}{d t} + \frac{P(t)}{R_2C} = \frac{Q(t)}{C} \left (1 + \frac{ R_1}{R_2} \right ) + R_1 \frac{d Q(t)}{dt} \quad \text{ (the misspesified model we use to fit the data) }\text{ [WK3] }\\ \eta: \quad & Q(t) = \frac{1}{R}P(t) + C \frac{dP(t)}{dt} \quad \text{ (the model we use to simulate data) } \text{ [WK2] } \end{align}\]

# uncomment to install
# install.packages("rstan")
# install.packages("ggplot2")
# install.packages("tidyverse")
library(rstan)
library(ggplot2)
library(tidyverse)
rstan_options(auto_write = TRUE)
options(mc.cores = 3) # allocate 3 cores (for each model we run 3 chains in parallel)
# Numerical simulator of the WK3 model
source("functions/WK2and3_sim_fn.R")
# Load flow data 
d = readRDS("Data/Inflow_time.rds")

Rtrue = 1; Ctrue = 1.1; Ztrue = 0.05 
flow = d$inflow*0.95
time = d$time
nP = 90 # number of pressure data
nI = 100# number of inflow data
nc = 1  # number of cardiac cycles
nflow = length(flow)
# 1. simulate WK3 data (R=R_2, Z=R_1)
Psim = WK3_simulate(flow = flow, time = time, R = Rtrue, C = Ctrue, Z=Ztrue) # simulate WK3 data for a given flow Q(t)
P_true = Psim
# 2. choose pressure and inflow indices
indP = round(seq(1, nflow, length.out = nP)); indI = round(seq(1, nflow, length.out = nI))
yP_real = Psim[indP]; yI_real = flow[indI] # noise free fimulated pressure and flow
# 3. Add noise
# set.seed(0)
set.seed(1)
Pnoise = rnorm(nP*nc, 0, 4) # sample pressure noise from N(0, 4^2)
Inoise = rnorm(nI*nc, 0, 10) # sample flow noise from N(0,10^2)
yP_real = rep(yP_real,nc) 
yI_real = rep(yI_real,nc)
# 4. store individual data in the population matrices
yP = yP_real + Pnoise # add noise
yI = yI_real + Inoise # add noise
tP = time[indP] # corresponding time (synchronized for the two cycles)
tI = time[indI] # corresponding time (synchronized for the two cycles)
data_PI = list(nP=nc*nP, nI=nc*nI, tP=rep(tP,nc), tI=rep(tI,nc), yP=yP, yI=yI, mP=12, mI=10)

WK2_VFE = stan_model('STAN/Approximations/VFE/WK2_delta_VFE.stan')
kp = kmeans(data.frame(x=data_PI$tP), centers = data_PI$mP)
ki = kmeans(data.frame(x=data_PI$tI), centers = data_PI$mI)
init = list("zP" = as.vector(kp$centers), "zI" = as.vector(ki$centers))
op_VFE=optimizing(WK2_VFE, data=data_PI, hessian=FALSE, init = init, verbose=TRUE, seed=0)

Chain 1: Initial log joint probability = -13400
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       19      -962.859     0.0138579        2.3978           1           1       31   
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       39      -949.581     0.0272867       2.43968           1           1       56   
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       49      -949.485   0.000188747     0.0088334           1           1       68   
Chain 1: Optimization terminated normally: 
Chain 1:   Convergence detected: relative gradient magnitude is below tolerance

op_VFE

$par
         rho        alpha        rho_d      alpha_d       mu_wk2       sigmaP       sigmaI            R            C        zP[1]        zP[2]        zP[3]        zP[4] 
 0.195718978  4.897341220  0.057918221  0.001593782 97.398040198  9.760864357 63.624564498  0.955552758  1.013489858  0.080072408  0.133867715  0.455149871  0.194511717 
       zP[5]        zP[6]        zP[7]        zP[8]        zP[9]       zP[10]       zP[11]       zP[12]        zI[1]        zI[2]        zI[3]        zI[4]        zI[5] 
 0.274922964  0.724878480  0.826984618  0.944283254  0.365017283  0.544772397  0.020043511  0.634880324  0.970363221  0.614334337  0.714632658  0.049797442  0.394214611 
       zI[6]        zI[7]        zI[8]        zI[9]       zI[10] 
 0.511302773  0.809613912  0.274381891  0.161115333  0.894987597 

$value
[1] -949.4846

$return_code
[1] 0

$theta_tilde
          rho    alpha      rho_d     alpha_d   mu_wk2   sigmaP   sigmaI         R       C      zP[1]     zP[2]     zP[3]     zP[4]    zP[5]     zP[6]     zP[7]     zP[8]
[1,] 0.195719 4.897341 0.05791822 0.001593782 97.39804 9.760864 63.62456 0.9555528 1.01349 0.08007241 0.1338677 0.4551499 0.1945117 0.274923 0.7248785 0.8269846 0.9442833
         zP[9]    zP[10]     zP[11]    zP[12]     zI[1]     zI[2]     zI[3]      zI[4]     zI[5]     zI[6]     zI[7]     zI[8]     zI[9]    zI[10]
[1,] 0.3650173 0.5447724 0.02004351 0.6348803 0.9703632 0.6143343 0.7146327 0.04979744 0.3942146 0.5113028 0.8096139 0.2743819 0.1611153 0.8949876

FITC

zP_opt_VFE=op_VFE$par[grep("zP",names(op_VFE$par))]
zI_opt_VFE=op_VFE$par[grep("zI",names(op_VFE$par))]
# plot(sort(zP_opt_VFE))
# plot(sort(zI_opt_VFE))
data_PI_Z_VFE= data_PI
data_PI_Z_VFE$zP = zP_opt_VFE
data_PI_Z_VFE$zI = zI_opt_VFE
fit_post_VFE=stan(file='STAN/Approximations/VFE/WK2_delta_VFE_fixed_Z.stan',
                  data=data_PI_Z_VFE,
                  chains=3,
                  iter=1000,
                  seed=0
)

starting worker pid=67458 on localhost:11354 at 11:30:29.975
starting worker pid=67472 on localhost:11354 at 11:30:30.172
starting worker pid=67486 on localhost:11354 at 11:30:30.365

SAMPLING FOR MODEL 'WK2_delta_VFE_fixed_Z' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.009392 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 93.92 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'WK2_delta_VFE_fixed_Z' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.010761 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 107.61 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 

SAMPLING FOR MODEL 'WK2_delta_VFE_fixed_Z' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.009607 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 96.07 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 58.2167 seconds (Warm-up)
Chain 3:                47.0122 seconds (Sampling)
Chain 3:                105.229 seconds (Total)
Chain 3: 
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 60.1515 seconds (Warm-up)
Chain 2:                47.3321 seconds (Sampling)
Chain 2:                107.484 seconds (Total)
Chain 2: 
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 69.9765 seconds (Warm-up)
Chain 1:                51.3671 seconds (Sampling)
Chain 1:                121.344 seconds (Total)
Chain 1: 

# stan_hist(fit_post_VFE)
stan_trace(fit_post_VFE)

zP_opt=op_FITC$par[grep("zP",names(op_FITC$par))]
zI_opt=op_FITC$par[grep("zI",names(op_FITC$par))]
data_PI_Z_FITC = data_PI
data_PI_Z_FITC$zP = zP_opt
data_PI_Z_FITC$zI = zI_opt
# plot(sort(zP_opt))
# plot(sort(zI_opt))

zP_opt=op_FITC$par[grep("zP",names(op_FITC$par))]
zI_opt=op_FITC$par[grep("zI",names(op_FITC$par))]
data_PI_Z_FITC = data_PI
data_PI_Z_FITC$zP = zP_opt
data_PI_Z_FITC$zI = zI_opt
# plot(sort(zP_opt))
# plot(sort(zI_opt))

fit_post_FITC=stan(file='STAN/Approximations/FITC/WK2_delta_FITC_fixed_Z.stan',
                  data=data_PI_Z_FITC,
                  chains=3,
                  iter=1000,
                  seed=0
)

starting worker pid=67510 on localhost:11354 at 11:32:59.213
starting worker pid=67524 on localhost:11354 at 11:32:59.413
starting worker pid=67538 on localhost:11354 at 11:32:59.614

SAMPLING FOR MODEL 'WK2_delta_FITC_fixed_Z' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.007932 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 79.32 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'WK2_delta_FITC_fixed_Z' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.007228 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 72.28 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'WK2_delta_FITC_fixed_Z' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.006141 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 61.41 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 60.4329 seconds (Warm-up)
Chain 1:                30.205 seconds (Sampling)
Chain 1:                90.6379 seconds (Total)
Chain 1: 
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 56.8 seconds (Warm-up)
Chain 3:                38.3528 seconds (Sampling)
Chain 3:                95.1527 seconds (Total)
Chain 3: 
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 62.991 seconds (Warm-up)
Chain 2:                38.5151 seconds (Sampling)
Chain 2:                101.506 seconds (Total)
Chain 2: 

# stan_hist(fit_post_FITC)
stan_trace(fit_post_FITC)

post_VFE = data.frame(rstan::extract(fit_post_VFE))
post_FITC = data.frame(rstan::extract(fit_post_FITC))
pr = c("R", "C", "sigmaP", "sigmaI")
pVFE = as.vector(as.matrix(post_VFE[,pr]))
pFITC = as.vector(as.matrix(post_FITC[,pr]))
Ns = nrow(post_VFE)
df_plot_post = data.frame(post = c(pVFE, pFITC), par = rep(rep(pr, each = Ns),2), model = c(rep("VFE",length(pr)*Ns), rep("FITC",length(pr)*Ns)))

mod_dat = df_plot_post%>%
  mutate(par = recode(par,
    "R" = "R",
    "C" = "C",
    "sigmaP" = "sigma[P]",
    "sigmaI" = "sigma[Q]"
  ))
df_true_par = data.frame(post=c(Rtrue+Ztrue, Ctrue, 4, 10),par = c("R", "C", "sigma[P]", "sigma[Q]"))
set_lim = data.frame(x=c(0.5,3.1),y=c(600, 600))
df_point_est = data.frame(val=c(op_VFE$par[pr],op_FITC$par[pr]),par = rep(df_true_par$par,2), model = rep(rep(c("VFE","FITC"),each=4),2))
pl_post = ggplot()+
  geom_histogram(data = mod_dat, aes(x=post, fill = par), color="black",bins = 60)+
  facet_grid(model~par, scales = "free", labeller = labeller(par = label_parsed))+
  geom_point(data = set_lim, aes(x=x,y=y), color = "white",alpha=0.000001, size=0.000001)+ # set limits on R and C
  geom_vline(aes(xintercept = post, linetype = "true"), data=df_true_par, size=0.8)+
  geom_vline(aes(xintercept = val, linetype = "map"), data=df_point_est, size=0.8)+
  theme_bw()+ 
  theme(#legend.position = "none",
        legend.title = element_blank(),
        legend.position="bottom",
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        strip.text.x = element_text(size = 13),
        strip.text.y = element_text(size = 13))+ 
  xlab("") + ylab("")+
  scale_fill_manual(
  breaks=c("R", "C", "sigma[P]", "sigma[Q]"),
  values=c("#F8766D","#00BE67","white", "white"),guide = "none")
(pl_post=pl_post + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))

ggsave("figures/Appr_delta_post.pdf", plot = pl_post, width = 20, height = 12, units = "cm")

nP_pred = 40
ind_P_pred = round(seq(1,length(time),length.out = nP_pred))
tP_pred = time[ind_P_pred]
data_PI_Z_FITC$nP_pred = nP_pred
data_PI_Z_FITC$tP_pred = tP_pred
data_PI_Z_FITC$nI_pred = nP_pred
data_PI_Z_FITC$tI_pred = tP_pred
N_samples = nrow(post_FITC)
data_post_FITC = list(alpha=post_FITC$alpha, rho=post_FITC$rho, alpha_d=post_FITC$alpha_d
                   , rho_d=post_FITC$rho_d, sigmaP=post_FITC$sigmaP, sigmaI=post_FITC$sigmaI
                   , R=post_FITC$R, C=post_FITC$C, N_samples=N_samples
  )
data_pred_FITC = c(data_PI_Z_FITC, data_post_FITC)

pred_FITC = stan(file = 'STAN/Approximations/FITC/FITC_delta_predictions.stan',
                 data = data_pred_FITC,
                 chains = 1, iter = 1, seed=123,
                 algorithm = "Fixed_param")

SAMPLING FOR MODEL 'FITC_delta_predictions' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0 seconds (Warm-up)
Chain 1:                3.00098 seconds (Sampling)
Chain 1:                3.00098 seconds (Total)
Chain 1: 

post_mu_CIs.fn = function(post_pred, ci = c(0.05,0.95), time = tP_pred){
  pp = rstan::extract(post_pred)
  dfP = pp$y_P[1,,]
  Psmr = data.frame(mean = colMeans(dfP), 
                    lower = apply(dfP, 2, quantile, probs = ci[1]), 
                    upper = apply(dfP, 2, quantile, probs = ci[2]),
                    time=time)
  dfI = pp$y_I[1,,]
  Ismr = data.frame(mean = colMeans(dfI), 
                    lower = apply(dfI, 2, quantile, probs = ci[1]), 
                    upper = apply(dfI, 2, quantile, probs = ci[2]),
                    time=time)
  return(list(Psmr=Psmr, Ismr=Ismr))
}

data_post_VFE = list(alpha=post_VFE$alpha, rho=post_VFE$rho, alpha_d=post_VFE$alpha_d
                      , rho_d=post_VFE$rho_d, sigmaP=post_VFE$sigmaP, sigmaI=post_VFE$sigmaI
                      , R=post_VFE$R, C=post_VFE$C, N_samples=N_samples
)
data_pred_VFE = c(data_PI_Z_VFE, data_post_VFE)
data_pred_VFE$nP_pred = nP_pred
data_pred_VFE$tP_pred = tP_pred
data_pred_VFE$nI_pred = nP_pred
data_pred_VFE$tI_pred = tP_pred

pred_VFE = stan(file = 'STAN/Approximations/VFE/WK2_delta_VFE_predictions.stan',
                 data = data_pred_VFE,
                 chains = 1, iter = 1, seed=123,
                 algorithm = "Fixed_param")

SAMPLING FOR MODEL 'WK2_delta_VFE_predictions' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0 seconds (Warm-up)
Chain 1:                2.47053 seconds (Sampling)
Chain 1:                2.47053 seconds (Total)
Chain 1: 

pp_VFE = rbind(post_mu_CIs.fn(post_pred=pred_VFE)$Psmr, post_mu_CIs.fn(post_pred=pred_VFE)$Ismr)
pp_VFE$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_VFE$model = "VFE"
pp_FITC = rbind(post_mu_CIs.fn(post_pred=pred_FITC)$Psmr,post_mu_CIs.fn(post_pred=pred_FITC)$Ismr)
pp_FITC$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_FITC$model = "FITC"
pred_df = rbind(pp_VFE, pp_FITC)
# head(pred_df)

Plug in noise estimates

We observe that the VFE model can severely overestimate the noise and therefore result in underfitting. A remedy to this problem is to fix the noise parameter of the functions \(P(t)\) and \(Q(t),\) (\(\sigma_P\) and \(\sigma_I\)). A possible solution for obtaining estimates for the noise parameters is to fit a standard GP model for each dataset \((y_P,t_P)\) and \(y_Q, t_Q\) indepentently and obtain MLE estimates via maximizing the marginal log-likelihood.

df_zP_VFE = data.frame(z=data_PI_Z_VFE$zP, y = rep(60, data_PI_Z_VFE$mP), model = "VFE", output = "pressure (mmHg)")
df_zI_VFE = data.frame(z=data_PI_Z_VFE$zI, y = rep(-100, data_PI_Z_VFE$mI), model = "VFE", output = "inflow (ml/min)")
df_zP_FITC = data.frame(z=data_PI_Z_FITC$zP, y = rep(60, data_PI_Z_FITC$mP), model = "FITC", output = "pressure (mmHg)")
df_zI_FITC = data.frame(z=data_PI_Z_FITC$zI, y = rep(-100, data_PI_Z_FITC$mI), model = "FITC", output = "inflow (ml/min)")
df_z = rbind(df_zP_VFE, df_zI_VFE, df_zP_FITC, df_zI_FITC)
P_true = data.frame(val=Psim, time=time)
P_true$output = "pressure (mmHg)"
I_true = data.frame(val=flow, time=time)
I_true$output = "inflow (ml/min)"
true_out = rbind(P_true, I_true)
obsP = data.frame(value=data_PI$yP, time = data_PI$tP, output = "pressure (mmHg)")
obsI = data.frame(value=data_PI$yI, time = data_PI$tI, output = "inflow (ml/min)")
obs = rbind(obsP,obsI)
pl_pred=ggplot()+
  geom_point(data = obs, aes(y=value, x=time, colour = "observed"), shape = 20)+
  geom_line(data = pred_df, aes(y=mean, x=time, linetype = "mean"), size=0.9)+
  geom_line(data = true_out, aes(y=val, x=time, linetype="true"), size=0.9)+
  geom_ribbon(data = pred_df,aes(ymin=lower, ymax=upper, x=time, fill = "90% CI"), alpha = 0.3)+
  facet_grid(output~model,scales = "free")+
  geom_point(data = df_z, aes(x=z, y=y, shape="Z"), size=3)+
  scale_fill_manual("",values=c("90% CI" = "grey12"))+
  theme_bw()+xlab("time (sec)")+ylab("")+
  scale_shape_manual("", values = c("Z" = 4))+
  theme(#legend.position = "none",
        legend.title = element_blank(),
        legend.position="bottom",
        strip.text.x = element_text(size = 13),
        strip.text.y = element_text(size = 10))
(pl_pred=pl_pred + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))

ggsave("figures/Appr_delta_pred.pdf", plot = pl_pred, width = 16, height = 10, units = "cm")

nsP = 25
indP = seq(1,data_PI$nP,length.out = nsP)
data_sample_P = list(N=nsP, x = data_PI$tP[indP], y = data_PI$yP[indP])
data_sample_I = list(N=nsP, x = data_PI$tI[indP], y = data_PI$yI[indP])
GP = stan_model('STAN/Approximations/GP_full/GP.stan')
GP_MLE_P=optimizing(GP, data=data_sample_P, hessian=FALSE, verbose=TRUE,seed=0)

Chain 1: Initial log joint probability = -84.3916
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       18       -65.767    0.00556334   7.66548e-05           1           1       30   
Chain 1: Optimization terminated normally: 
Chain 1:   Convergence detected: relative gradient magnitude is below tolerance

GP_MLE_P

$par
       rho      alpha      sigma 
 0.2147435 76.0190356  3.3560127 

$value
[1] -65.767

$return_code
[1] 0

$theta_tilde
           rho    alpha    sigma
[1,] 0.2147435 76.01904 3.356013

sigma_P_MLE = GP_MLE_P$par["sigma"] 
GP_MLE_I=optimizing(GP, data=data_sample_I, hessian=FALSE, verbose=TRUE,seed=0)

Chain 1: Initial log joint probability = -656.901
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       19      -111.136      0.558199      0.622198           1           1       26   
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       39      -111.126      0.321616   0.000908148      0.4931     0.04931       55   
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       41      -111.126      0.373032    0.00159772      0.1782           1       58   
Chain 1: Optimization terminated normally: 
Chain 1:   Convergence detected: relative gradient magnitude is below tolerance

GP_MLE_I

$par
         rho        alpha        sigma 
  0.07883712  99.99999925  11.38052032 

$value
[1] -111.1259

$return_code
[1] 0

$theta_tilde
            rho alpha    sigma
[1,] 0.07883712   100 11.38052

sigma_I_MLE = GP_MLE_I$par["sigma"] 

data_pred_VFE_MLE = data_pred_VFE
data_pred_VFE_MLE$sigmaP = sigma_P_MLE
data_pred_VFE_MLE$sigmaI = sigma_I_MLE
pred_VFE_MLE = stan(file = 'STAN/Approximations/VFE/MLE_sigma/WK2_delta_VFE_predictions.stan',
                 data = data_pred_VFE_MLE ,
                 chains = 1, iter = 1, seed=123,
                 algorithm = "Fixed_param")

SAMPLING FOR MODEL 'WK2_delta_VFE_predictions' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0 seconds (Warm-up)
Chain 1:                2.4886 seconds (Sampling)
Chain 1:                2.4886 seconds (Total)
Chain 1: 

pp_VFE_MLE = rbind(post_mu_CIs.fn(post_pred=pred_VFE_MLE)$Psmr, post_mu_CIs.fn(post_pred=pred_VFE_MLE)$Ismr)
pp_VFE_MLE$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_VFE_MLE$model = "VFE fixed noise"
pp_VFE = rbind(post_mu_CIs.fn(post_pred=pred_VFE)$Psmr, post_mu_CIs.fn(post_pred=pred_VFE)$Ismr)
pp_VFE$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_VFE$model = "VFE"
pred_df = rbind(pp_VFE_MLE,pp_VFE)
df_zP_VFE = data.frame(z=data_PI_Z_VFE$zP, y = rep(60, data_PI_Z_VFE$mP), model = "VFE", output = "pressure (mmHg)")
df_zI_VFE = data.frame(z=data_PI_Z_VFE$zI, y = rep(-100, data_PI_Z_VFE$mI), model = "VFE", output = "inflow (ml/min)")
df_zP_VFE_MLE = data.frame(z=data_pred_VFE_MLE$zP, y = rep(60, data_pred_VFE_MLE$mP), model = "VFE fixed noise", output = "pressure (mmHg)")
df_zI_VFE_MLE = data.frame(z=data_pred_VFE_MLE$zI, y = rep(-100, data_pred_VFE_MLE$mI), model = "VFE fixed noise", output = "inflow (ml/min)")
df_z = rbind(df_zP_VFE, df_zI_VFE,df_zP_VFE_MLE,df_zI_VFE_MLE)
P_true = data.frame(val=Psim, time=time)
P_true$output = "pressure (mmHg)"
I_true = data.frame(val=flow, time=time)
I_true$output = "inflow (ml/min)"
true_out = rbind(P_true, I_true)
obsP = data.frame(value=data_PI$yP, time = data_PI$tP, output = "pressure (mmHg)")
obsI = data.frame(value=data_PI$yI, time = data_PI$tI, output = "inflow (ml/min)")
obs = rbind(obsP,obsI)
pl_pred=ggplot()+
  geom_point(data = obs, aes(y=value, x=time, colour = "observed"), shape = 20)+
  geom_line(data = pred_df, aes(y=mean, x=time, linetype = "mean"), size=0.9)+
  geom_line(data = true_out, aes(y=val, x=time, linetype="true"), size=0.9)+
  geom_ribbon(data = pred_df,aes(ymin=lower, ymax=upper, x=time, fill = "90% CI"), alpha = 0.3)+
  facet_grid(output~model,scales = "free")+
  geom_point(data = df_z, aes(x=z, y=y, shape="Z"), size=3)+
  scale_fill_manual("",values=c("90% CI" = "grey12"))+
  theme_bw()+xlab("time (sec)")+ylab("")+
  scale_shape_manual("", values = c("Z" = 4))+
  theme(#legend.position = "none",
        legend.title = element_blank(),
        legend.position="bottom",
        strip.text.x = element_text(size = 13),
        strip.text.y = element_text(size = 10))
(pl_pred=pl_pred + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))

ggsave("figures/Appr_delta_pred_noise.pdf", plot = pl_pred, width = 16, height = 10, units = "cm")